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Abstract 

This  article  is  the  second  of  a  two-part  paper  on  ELCIRC,  an  Eulerian-Lagrangian  finite  difference/finite  volume  model 
designed  to  simulate  3D  baroclinic  circulation  across  river-to-ocean  scales.  In  part  one  (Zhang  et  al.,  2004),  we  described 
the  formulation  of  ELCIRC  and  assessed  its  baseline  numerical  skill.  Here,  we  describe  the  application  of  ELCIRC  within 
CORIE,  a  coastal  margin  observatory  for  the  Columbia  River  estuary  and  plume.  We  first  introduce  the  CORIE  modeling 
system  and  its  multiple  modes  of  simulation,  external  forcings,  observational  controls,  and  automated  products.  We  then 
focus  on  the  evaluation  of  highly  resolved,  year-long  ELCIRC  simulations,  using  two  variables  (water  level  and  salinity)  to 
illustrate  simulation  quality  and  sensitivity  to  modeling  choices.  We  show  that,  process-wise,  simulations  capture  well 
important  aspects  of  the  response  of  estuarine  and  plume  circulation  to  ocean,  river,  and  atmospheric  forcings. 
Quantitatively,  water  levels  are  robustly  represented,  while  salinity  intrusion  and  plume  dynamics  remain  challenging.  Our 
analysis  highlights  the  benefits  of  conducting  model  evaluations  over  large  time  windows  (months  to  years),  to  avoid 
significant  localized  biases.  The  robustness  and  computational  efficiency  of  ELCIRC  has  proved  invaluable  in  identifying 
and  reducing  non-algorithmic  sources  of  errors,  including  parameterization  (e.g.,  turbulence  closure  and  stresses  at  the  air- 
water  interface)  and  external  forcings  (e.g.,  ocean  conditions  and  atmospheric  forcings). 
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1.  Introduction 

Integrated  observatories  (Clark  and  Isern,  2003; 
Martin,  2003;  USCOP,  2004)  are  expected  to 
dramatically  improve  the  understanding  of  the 
ocean  across  scales  and  processes,  and  to  provide 
unprecedented,  objective  information  to  address 
societal  priorities  regarding  ocean  preservation 
and  utilization.  Meeting  these  lofty  expectations 
will  require  improvements  in  underlying  technol¬ 
ogies  (e.g.,  models,  platforms,  sensors,  and  in¬ 
formation  technologies),  as  well  as  adjustments  in 
their  use. 

A  preview  of  challenges  to  come  has  been 
provided  by  selected  prototype  ocean  observa¬ 
tories  (Parker,  1997;  Glenn  et  al.,  2000;  Steere 
et  al.,  2000;  Rhodes  et  al.,  2001;  Baptista,  2002). 
Developed  and  maintained  since  1996  for  the 
Columbia  River  estuary  and  plume,  CORIE 
(CCALMR,  1996-2004;  Baptista  et  al.,  1998; 
Baptista  et  al.,  1999)  is  one  such  prototype. 
CORIE  was  designed  from  the  onset  as  a  multi¬ 
purpose,  regional  infrastructure  for  research, 
education,  and  management.  The  design  includes 
three  integrated  components:  a  real-time  observa¬ 
tion  network,  a  modeling  system,  and  a  web-based 
information  system. 

Perhaps  surprisingly,  of  the  three  CORIE 
components,  the  modeling  system  has  posed  the 
most  fundamental  challenges,  calling  for  new 
modeling  technologies  and  paradigms.  In  particu¬ 
lar,  we  found  the  need  to  develop  a  new  cross-scale 
3D  baroclinic  circulation  code  (ELCIRC;  Zhang 
et  al.,  2004)  in  order  to  meet  operational  require¬ 
ments  of  efficiency  and  quality.  Also,  automated 
integrative  procedures — including  the  generation 
of  model  forcings,  quality  controls  and  modeling 
products — have  become  essential  to  creating, 
improving,  and  interpreting  simulations.  More¬ 
over,  multiple  simulation  modes — including  daily 
forecasts,  multi-year  hindcasts,  and  scenario  simu¬ 
lations — forced  the  development  of  multi-scale, 
long-term  calibration  and  validation  strategies  and 
procedures. 

Here,  we  describe  the  CORIE  modeling  system 
(Section  2)  and  present  selected  results  (Section  3), 
with  emphasis  on  water  levels  and  salinities.  The 
description  of  the  modeling  system  is  intended  as  a 


reference  for  derivative  papers,  which  will  comple¬ 
ment  the  present  work  by  exploring  in  further 
depth  specific  modeling  aspects,  and  scientific  and 
management  applications  of  CORIE.  The  results 
in  Section  3  illustrate  the  extent  to  which  the 
modeling  system  is  already  able  to  describe 
complex,  multi-scale  circulation  processes  in  the 
Columbia  River  and  help  identify  directions  for 
further  improvement.  In  Section  4,  we  discuss 
implications  for  derivative  research  on  Columbia 
River  circulation,  for  further  algorithmic  develop¬ 
ment  of  cross-scale  numerical  models  and  for 
coastal  estuarine  and  plume  modeling  within  the 
requirements  of  ocean  observatories. 


2.  The  CORIE  modeling  system 

One  of  the  world’s  classic  river-dominated 
estuaries,  the  Columbia  River,  is  a  highly  dynamic 
system  that  responds  dramatically  to  changes  in 
ocean  tides  and  water  properties,  regulated  river 
discharges,  and  coastal  winds.  A  dominant  hydro- 
graphic  feature  on  the  U.S.  west  coast,  the 
Columbia  River  plume  exports  dissolved  and 
particulate  matter  hundreds  of  kilometers  along 
and  across  the  continental  shelf  (Barnes  et  al., 
1972;  Grimes  and  Kingsford,  1996).  In  response  to 
seasonal  changes  of  the  large-scale  coastal  circula¬ 
tion  patterns,  the  plume  typically  develops  north¬ 
ward  along  the  coastal  shelf  in  fall  and  winter,  and 
southwestward  offshore  of  the  shelf  in  spring  and 
summer.  However,  the  direction,  thickness,  and 
width  of  the  plume — and,  in  particular,  the  near¬ 
field  plume — can  change  in  hours  to  days  in 
response  to  local  winds  (CCALMR,  1996-2004; 
Hickey  et  al.,  1998;  Garcia-Berdeal  et  al.,  2002). 

Compressed  and  often  stratified,  the  Columbia 
River  estuary  is  subject  to  extreme  variations  in 
salinity  intrusion  and  stratification  regime  (Ha¬ 
milton,  1990;  Jay  and  Smith,  1990;  Chawla  et  al., 
in  prep.).  Two  main  channels  (one  dredged 
for  navigation)  cut  the  otherwise  shallow  estuary 
(Fig.  1).  A  shallow  coastal  region  north  of  the 
Columbia  River  mouth  combines  with  Coriolis  to 
establish  an  underlying  tendency  for  the  near¬ 
plume  to  move  north.  This  northward  tendency  is 
countered  by  the  exit  angle  of  the  navigation 
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Fig.  1.  The  North  Channel  and  Navigation  Channel  cut  the  otherwise  shallow  Columbia  River  estuary.  Lateral  bays  are  shallow  and 
ecologically  important  environments.  Shallow  bathymetry  immediately  North  of  the  mouth  of  the  estuary  combines  with  Coriolis  to 
naturally  bend  the  plume  northward  in  the  absence  of  winds. 


channel  and  is  either  countered  or  reinforced  by 
local  winds. 

The  CORIE  modeling  system  integrates  models, 
forcings,  and  quality  controls  to  produce  daily 
forecasts  and  multi-year  simulation  databases  of 
circulation  in  both  the  estuary  and  the  plume 
(Fig.  2).  Data  assimilation  is  used  only  sparingly 
(Section  2.3).  Elements  of  complexity  for  the 
simulations  include:  (a)  highly  variable  and  non¬ 
linear  forcings  (e.g.,  wind,  tides,  river  discharge); 
(b)  dynamic  density  fronts  and  strong  buoyancy- 
driven  flow;  and  (c)  broad  range  of  inter-connected 
spatial  (from  under  100  m  to  well  above  10  km) 
and  temporal  scales  (minutes  to  decades). 

To  address  these  complexities,  we  found  it 
necessary  to  develop  a  new  3D  baroclinic  circula¬ 
tion  model,  ELCIRC.  The  motivation,  formula¬ 
tion,  and  basic  skill  assessment  of  ELCIRC  have 
been  presented  in  Zhang  et  al.  (2004).  In  summary, 
ELCIRC  is  a  finite-difference/finite-volume  model, 
based  on  unstructured  grids;  the  numerical  scheme 
is  robust,  volume  conservative,  and  computation¬ 
ally  efficient.  Within  the  constraints  of  the  hydro¬ 
static  approximation,  the  ELCIRC  physical 
formulation  allows  for  a  unified  river-to-ocean 
approach,  allying  to  features  of  state-of-the-art 
coastal  ocean  models  the  ability  to  treat  estuary 


features  such  as  regions  of  high  advection  and  of 
extensive  wetting  and  drying. 

Taking  advantage  of  the  characteristics  of 
ELCIRC,  the  domain  of  simulation  spans  from 
the  Bonneville  Dam  and  the  Willamette  Falls, 
approximately  240  km  upstream  of  the  entrance  to 
the  Columbia  River  estuary,  to  and  beyond  the 
continental  shelf  of  California,  Oregon,  Washing¬ 
ton,  and  British  Columbia.  Included  in  the 
domain,  albeit  at  low  resolution,  are  the  straits 
of  Juan  de  Fuca  and  Georgia,  and  the  Puget 
Sound.  Computational  grids  (Fig.  3;  see  also 
Figs.  6a  and  15c  for  local  details)  are  typically 
3D,  unstructured  in  the  horizontal  plane,  with  z- 
coordinates  in  the  vertical.  Computational  and 
archival  requirements  are  met  with  a  dedicated 
computer  infrastructure.  At  the  time  of  this 
writing,  the  infrastructure  includes  20  dual  CPU 
Intel  compute  nodes  (2.4  GHz,  4  Gb)  organized  as 
a  Beowulf  cluster  and  28  TB  of  online  disk  arrays. 
For  a  baroclinic  simulation  on  a  typical  computa¬ 
tional  grid  (34190  horizontal  nodes;  50622  hor¬ 
izontal  hybrid  elements;  62  z-levels)  and  with  a 
typical  time  step  (1.5  minutes),  ELCIRC  runs 
2.5-3  times  faster  than  real  time  in  a  single  CPU. 
The  upper  bound  in  this  range  applies  to  simula¬ 
tions  with  a  zero-equation  turbulence  closure,  and 


938 


A.M.  Baptistci  et  al.  /  Continental  Shelf  Research  25  (2005)  935-972 


Fig.  2.  The  CORIE  modeling  system  integrates  models  with  external  forcings,  quality  controls,  and  products. 


the  lower  bound  applies  to  simulations  with  a 
2.5-equation  turbulence  closure.  Purely  barotropic 
simulations  can  be  run  with  a  1 5-min  time  step  and 
are  about  60  times  faster  than  real-time.  However, 
rarely  is  it  justified  to  run  a  barotropic  simulation 
in  the  Columbia  River,  due  to  strong  buoyancy 
effects.  About  0.8  TB  of  online  storage  is  required 
to  store  a  typical  year-long  baroclinic  simulation. 

2.1.  External  forcings 

Ocean,  atmospheric,  and  river  influences  control 
much  of  the  dynamics  of  the  Columbia  River 
estuary  and  plume.  All  of  these  influences  are 
highly  variable  in  space  and  time.  It  is  essential 
to  account  for  this  variability,  for  a  successful 
simulation  of  circulation.  We  describe  below 
strategies  and  information  sources  to  characterize 
external  forcings  in  the  CORIE  modeling 
system.  Where  they  diverge,  we  briefly  note 


differences  between  strategies  to  address  forcings 
for  daily  forecasts  and  strategies  for  retrospective 
simulations  (hindcast  databases  and  calibration 
runs). 

2.1.1.  River  inputs 

The  Columbia  River  and  the  Fraser  River 
watersheds  are  the  dominant  freshwater  source 
for  the  Pacific  Northwest  (see  Fig.  3  for  geographic 
reference),  and  both  are  accounted  for  in  the 
CORIE  modeling  system.  Within  the  Columbia 
River,  freshwater  inputs  are  considered  from  the 
Bonneville  Dam  (for  the  main  stem  of 
the  Columbia  River)  and  from  Newberg  (for  the 
Willamette  River).  Neglected  are  smaller  fresh¬ 
water  inputs,  including  those  from  the  Cowlitz, 
Lewis,  and  Sandy  rivers.  Bonneville  discharges 
retain — despite  heavy  hydropower  regulation — 
substantial  seasonal  (e.g.,  spring  freshets)  and 
inter-annual  variability  (e.g.,  in  response  to  El 
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Strait  of 


Fig.  3.  The  computational  domain  extends  from  California,  through  Oregon  and  Washington,  to  British  Columbia.  Grids  are  hybrid 
and  resolution  is  finest  in  the  estuary  and  near  plume  (e.g.,  see  details  in  Figs.  6a  and  1 5c). 

Nino— Southern  Oscillation  and  in  response  to  (US  Geological  Survey)  gauges:  USGS-ID 

Pacific — Decadal  Oscillation).  14128870,  downstream  of  the  Bonneville  dam; 

For  retrospective  simulations,  information  on  and  USGS-ID  14197900,  at  Newberg.  Tempera- 

freshwater  discharge  is  obtained  from  two  USGS  ture  data  for  Columbia  River  are  obtained  from 
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the  same  site  that  records  the  river  discharge. 
Temperature  data  for  Willamette  River  are  ob¬ 
tained  from  a  USGS  gauge  in  Portland  (USGS-ID 
14211720),  located  about  60  km  downstream  of 
the  Newberg  gauge.  Data  from  the  same  gauges 
are  used  for  the  daily  forecasts  through  short-term 
extrapolation.  Fraser  discharge  is  taken  from  the 
Hope  gauge  (08MF005)  of  Environment  Canada. 
For  forecasts,  the  actual  discharge  values  are 
unavailable,  so  they  are  replaced  with  a  discharge 
climatology  based  on  1912-2003  data.  Fraser 
temperature  is  extracted  from  the  NCOM  model 
results  (see  below). 

2.1.2.  Ocean  conditions 

We  have  used  four  semi-diurnal  (M2,  S2,  N2, 
K2)  and  four  diurnal  (Kj,  Oj,  Pi,  Ch)  constituents 
to  characterize  ocean  tides  at  the  boundaries  of  the 
CORIE  computational  domain.  Spatially  variable 
tidal  amplitudes  and  phases  are  available  for  the 
Eastern  North  Pacific  from  at  Egbert  et  al.  (1994) 
and  Myers  and  Baptista  (2001).  We  typically  use 
Myers  and  Baptista  (2001)  (e.g.,  Table  1)  after 
early  tests  showed  limited  sensitivity  to  the  choice 
of  one  or  another  source.  Nodal  factors  and 
astronomical  arguments  are  calculated  from  the 
tidal  package  of  Foreman  (1977).  No  low- 
frequency  set-up  was  imposed  at  the  boundaries. 
The  earth  tidal  potential  (from  Reid,  1990)  can  be 
included  but  is  typically  neglected. 

Seasonal  climatology  for  ocean  salinity  and 
temperature  is  available  from  Levitus  (1982). 
However,  Levitus  climatology  is  notably  coarse 
in  particular  for  coastal  applications  and  is 
inherently  unable  to  capture  the  inter-annual 
variability  of  salinity  and  temperature.  As  an 
alternative  to  climatology,  we  use  the  output  from 
operational  Navy  products  to  define  ocean  salinity 
and  temperature  conditions,  using  2002  hindcast 
database  simulations  as  the  benchmark.  Devel¬ 
oped  by  the  Naval  Research  Laboratory  for 
application  to  coastal  and  global  prediction  of 
ocean  dynamic  and  thermodynamic  fields  (Martin, 
2000),  the  Navy  Coastal  Ocean  Model  (NCOM)  is 
a  variant  of  the  Princeton  Ocean  Model  (POM; 
Blumberg  and  Mellor,  1987).  NCOM  runs  using 
1  /8th  degree  horizontal  resolution  and  40  vertical 
levels  that  are  a  combination  of  sigma  levels  in  the 


upper  150  m  of  the  ocean  and  z-levels  from  150  m 
to  the  ocean  bottom  (Rhodes  et  al.,  2001).  This 
model  has  been  tested  using  both  climatological 
forcings  and  real  time  atmospheric  forcings  from 
the  Navy’s  Operational  Global  Atmospheric  Pre¬ 
diction  System  (NOGAPS;  Rosmond  et  al.,  2002). 
Global  NCOM  has  been  spun  up  from  a 
climatological  state  to  the  present,  using  a 
combination  of  NOGAPS  forcings  and  the  assim¬ 
ilation  of  satellite  altimeter  data  and  3D  tempera¬ 
ture  and  salinity  observations  derived  from  the 
Modular  Ocean  Data  Assimilation  System 
(MODAS;  Fox  et  al.,  2002).  By  default,  we  use 
NCOM  results  to  define  the  initial  oceanic  salinity 
and  temperature  ( S ,  T)  conditions  as  well  as  for 
nudging  ELCIRC  S,  T  results  in  the  ocean 
(Section  2.3). 

To  date,  CORIE  simulations  have  relied  on  the 
internal  physics  of  the  ELCIRC  model  and  on 
external  (winds  and  atmospheric  pressure)  and 
internal  (density  gradients)  forcings  to  generate 
ocean  set-up  and  circulation  within  the  modeling 
domain. 

2.1.3.  Atmospheric  forcings 

Several  different  atmospheric  forcings  are  uti¬ 
lized  by  the  CORIE  modeling  system.  They  are 
divided  into  two  main  groups:  near-surface  atmo¬ 
spheric  properties  and  downwelling  radiative 
fluxes  at  the  surface.  The  atmospheric  properties 
include  the  x-  and  /-components  of  the  wind  at  a 
height  of  10  m,  the  surface  atmospheric  pressure 
(reduced  to  mean  sea  level),  the  air  temperature  at 
2  m,  and  the  specific  humidity  at  2  m.  The  radiative 
fluxes  include  the  downwelling  shortwave  (solar) 
and  the  downwelling  longwave  (infrared)  radia¬ 
tions  at  the  water  surface. 

These  forcing  data  have  been  compiled  from  a 
number  of  different  sources,  but  the  data  fall  into 
two  broad  divisions:  locally  archived  forecast  data 
and  re-analysis  data  (essentially  hindcast  simula¬ 
tions  with  extensive  data  assimilation  of  available 
observations).  The  forecast  datasets  include  data 
from:  (a)  the  National  Oceanic  and  Atmospheric 
Administration  (NOAA)  National  Centers  for 
Environmental  Prediction  (NCEP)  Global  Fore¬ 
cast  System  (GFS,  also  referred  to  as  the  MRF 
and/or  AVN  forecasts);  (b)  NCEP  Eta  Model;  and 
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(c)  the  Advanced  Regional  Prediction  System 
developed  at  the  University  of  Oklahoma,  as 
modified  and  run  at  Oregon  State  University 
(OSU/ARPS).  The  reanalysis  dataset  originates 
from  a  joint  project  of  NCEP  and  the  National 
Center  for  Atmospheric  Research  (NCAR;  the 
NCAR/NCEP  Global  Reanalysis  Project). 

2.2.  Observational  controls 

Inherent  to  the  concept  of  the  CORIE  modeling 
system  are  systematic  quality  controls  based  on 
comparisons  of  data  from  various  long-term 
observational  networks  (Figs.  4  and  5).  Most  of 
these  networks  have  real-time  or  quasi  real-time 
telemetry,  thus  supporting  both  hindcast  simula¬ 
tions  (databases  and  calibration  runs)  and  daily 
forecasts. 


The  CORIE  real-time  observation  network 
covers  the  estuary  extensively  and  the  near-ocean 
sparingly  (Fig.  4;  CCALMR  (1996-2004)).  At 
each  station,  in  situ  sensors  measure  various 
combinations  of  water  temperature,  salinity,  pres¬ 
sure,  velocity,  and  atmospheric  parameters.  Ob¬ 
servations  are  available  in  both  real-time  and  long¬ 
term  archives,  some  extending  nearly  8  years.  We 
have  also  equipped  the  M/V  Forerunner,  a  50-foot 
vessel,  with  a  pump-through,  conductivity-tem¬ 
perature  sensor  and  a  hull-mounted  acoustic 
Doppler  profiler.  The  vessel  is  also  used  for 
targeted  CORIE  cruises  (e.g.,  Fig.  20). 

Other  observation  networks  used  include 
NOAA  Center  for  Operational  Oceanographic 
Products  and  Services  (CO-OPS),  NOAA 
National  Data  Buoy  Center  (NDBC),  U.S. 
Army  Corps  of  Engineers  (USACE),  and  USGS 
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Fig.  4.  Fixed  stations  of  the  CORIE  observation  network  are  concentrated  on  the  estuary  up  to  the  limit  of  salinity  intrusion,  with  an 
offshore  presence  through  station  ogiOl.  Most  stations  are  in  piles  or  similar  structures,  but  Acoustic  Doppler  profilers  (available  at 
five  stations:  ogiO\,  red26,  tansy,  aml69,  and  am0\2)  require  either  bottom  (in  the  estuary)  or  surface  (shelf)  frames  or  buoys.  Besides 
fixed  stations,  CORIE  includes  a  mobile  station  in  the  form  of  a  training  vessel.  The  configuration  of  the  CORIE  network  and  its 
sensors  is  still  evolving;  see  (CCALMR,  1996-2004)  for  a  current  configuration.  The  planned  deployment  of  a  short  range  high- 
frequency  radar  and  an  x-band  radar  is  through  collaborations  with  researchers  Oregon  State  University. 
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Fig.  5.  Other  observation  networks  in  the  region.  NOAA  COOPS  stations  provide  tidal  elevation  data,  NDBC  buoys  provide 
atmospheric  winds  and  pressure  information,  and  USGS  and  USACE  gauges  provide  river  discharge  and/or  temperature  information. 


(cf.  Fig.  5).  Besides  systematic  quality  control  of 
CORIE  simulations  against  information  from 
long-term  fixed  observation  networks,  the  CORIE 
modeling  system  also  conducts  more  “random” 
quality  controls  against  episodic  data  from  a 
variety  of  sources.  Particularly  useful  are  Synthetic 
Aperture  Radar  satellite  images  (e.g.,  Fig.  15a) 
and  diverse  oceanographic  cruises. 

Key  oceanographic  cruises  include:  (a)  an 
extensive  1990-1991  plume  study  with  moorings 
and  vessels  (Hickey  et  al.,  1998);  (b)  extensive 
estuarine  cruises  in  the  Columbia  River  estuary 


conducted  as  a  part  of  the  Columbia  River 
component  (CRETM,  1990-2000)  of  a  national 
land-margin  ecosystem  research  program  (LMER 
Coordinating  Committee,  1992);  (c)  periodic 
CORIE  cruises  in  the  Columbia  River  estuary 
and  plume  (since  1 996);  and  (d)  cruises  conducted 
by  NOAA  Fisheries  since  1998  along  coastal 
transects  in  California,  Oregon,  and  Washington, 
and  at  or  near  the  mouth  of  the  Columbia  River. 
Recently,  extensive  Columbia  River  plume  surveys 
were  conducted  in  May-July  of  2004,  by  two 
different  but  overlapping  multi-investigator  teams, 
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using  diverse  and  sophisticated  observation  tech¬ 
niques  from  land,  airborne,  and  in-water  (moored 
and  mobile)  platforms  (Section  3.1). 

2.3.  Data  assimilation  and  nudging 

The  CORIE  philosophy  is  that  the  numerical 
circulation  model  is  ultimately  responsible  for 
representing  the  physics  inside  the  domain,  once 
appropriate  choices  have  been  made  for  para¬ 
meters  and  external  forcings.  Thus,  the  CORIE 
modeling  system  has  used  data  assimilation  spar¬ 
ingly:  primarily,  for  improving  the  definition  of 
external  forcings  and  for  off-line  optimization  of 
empirical  parameters. 

Specifically,  data  assimilation  was  used  to  define 
ocean  tidal  forcings  (Myers  and  Baptista,  2001) 
and  is  being  used  to  optimize  bottom  friction 
(Frolov  et  al.,  2004).  In  addition,  under  certain 
circumstances,  ocean  conditions  ( S ,  T)  are  locally 
nudged  to  information  from  either  global  circula¬ 
tion  models  or  climatology— which  in  turn,  have 
been  either  data-assimilated  or  objectively  ana¬ 
lyzed  from  data.  With  this  nudging,  we  seek  to 
impose  non-refiective  ocean  boundary  conditions 
and  to  moderate  errors  in  heat-balance  calcula¬ 
tions  resulting  from  specifying  imprecise  atmo¬ 
spheric  forcings. 

Nudging  is  accomplished  through  a  simple 
algorithm: 

Pn(x,y,z)  =  (1  -  a )FlELCIRC(x,y,z) 

+  aPbase(X’y’Z )  0) 

with 

a(x,y,  z)  =  y{x,y)\j/{z), 

where  /?"  is  the  nudged  value  at  time  n,  l>nELctRC 
the  value  computed  directly  by  solving  the 
governing  equations,  and  [Snbase  is  the  reference 
value  from  a  global  circulation  model  or  from 
climatology.  The  nudging  factor  y  is  typically 
zero  in  the  estuary  and  near-plume  but 
increases  toward  the  ocean  in  patterns  dictated 
by  the  objectives  of  the  particular  simulation.  The 
vertical  nudging  profile  ip  is  linear  or  piece-wise 
linear. 


2. 4.  Products 

The  large  number  of  simulations  conducted 
daily  within  CORIE  requires  systematic,  standar¬ 
dized  processing.  Automated  procedures  have 
been  established  for  separate  but  similar  proces¬ 
sing  of  daily  forecasts  and  of  hindcast  simulations 
(both  hindcast  databases  and  calibration  runs). 
The  results  of  hindcast  simulations  are  publicly 
accessible  (CCALMR,  1996-2004). 


3.  Selected  results 

Simulations  of  Columbia  River  circulation  are 
sensitive — often  in  complex,  nonlinear  manners — 
to  a  wide  range  of  modeling  choices,  including 
initialization  strategies  and  representation  of 
internal  parameters,  bathymetry,  and  external 
forcings. 

To  understand  and  to  address  this  sensitivity,  we 
have  developed  a  sustained,  iterative  strategy  that 
involves  simulations  at  multiple  time  windows. 
Anchoring  this  process  are  hindcast  circulation 
databases  that  extend  for  one  or  more  years  and 
that  take  2-3  months  to  generate.  Databases  reflect 
best-available  modeling  choices  at  the  beginning  of 
their  generation,  and  those  choices  are  kept 
unchanged  throughout  the  generation  of  the  entire 
database. 

Complementary  calibration  runs  are,  however, 
also  used  in  parallel  to  advance  the  state-of-the-art 
in  CORIE  simulations.  Thus,  it  is  rare  that  a 
database  is  completed  without  one  or  more 
modeling  choices  being  shown  to  be  non-opti- 
mal — and  this  information  feeds  the  next  data¬ 
base.  The  duration  of  calibration  runs  is  typically 
between  one  week  and  several  months. 

As  an  illustration  of  our  strategy,  in  Section  3.1 
we  will  describe  a  baseline  simulation  and  will  then 
summarize  how  alternative  modeling  options 
affect  the  representation  of  processes  and  vari¬ 
ables.  We  will  focus  our  analysis  on  just  two 
variables:  water  levels  and  salinity.  Velocities 
and  temperatures  will  be  addressed  in  separate 
publications.  However,  because  salt  propagation  is 
essentially  a  transport  problem,  the  present 
analysis  already  provides  a  stringent,  albeit 
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indirect,  assessment  of  the  circulation  capabilities 
of  the  numerical  model. 

3. 1.  Baseline  simulation 

We  take  as  a  reference  a  specific  circulation 
database  11  (DB11)  and  focus  on  its  results  for 
2002.  At  the  time  of  this  writing,  DB1 1  is  the  most 
current  and  comprehensive  CORIE  database, 
covering  1999-present.  As  typical  in  CORIE, 
DB11  modeling  choices  incorporate  lessons 
learned  in  several  previous  databases  and 
calibration  runs  (e.g.,  Section  3.2).  Shortcomings 
detected  during  the  generation  of  DB11  have 
already  inspired  parallel  calibrations  runs  (e.g.. 
Section  3.2).  Table  5  summarizes  parameters  used 
in  DB11,  as  well  as  in  other  databases  and 
calibration  runs  presented  in  this  paper. 

3.1.1.  Simulation  set-up 

DB11  was  constructed  by  combining  multiple 
simulation  ensembles  conducted  in  parallel.  Each 
ensemble  is  composed  of  about  3  months,  run 
sequentially.  To  minimize  discontinuities  across 
ensembles,  contiguous  ensembles  overlap  by  two 
(or  more  if  needed)  weeks:  the  first  two  weeks  of 
each  ensemble  are  considered  warm-up  and  are 
eventually  discarded  in  favor  of  the  last  weeks  of 
the  previous  ensemble.  The  rationale  behind  this 
“temporal”  parallelization  hinges  on  the  hypoth¬ 
esis  that  a  dynamic  equilibrium  is  established 
within  each  ensemble.  This  hypothesis  will  be 
tested  below.  All  simulations  include  nudging  of 
ocean  S,  T  to  NCOM  results,  which  also  provide 
initial  conditions  for  each  simulation  ensemble. 

The  computational  grid  uses  a  combination  of 
quadrangles  and  triangles,  with  highest  spatial 
resolution  concentrated  in  the  Columbia  River 
estuary  and  near  plume.  Over  35%  of  the  elements 
in  the  grid  have  an  equivalent  diameter  between 
100  and  200  m  (Fig.  6).  Upstream  of  the  estuary, 
the  main  river  channel  is  carefully  delineated,  but 
flood  plains  are  often  under-detailed.  Thus,  that 
part  of  the  grid  functions  primarily  to  delineate  a 
conduit  for  freshwater  discharge  into  the  estuary, 
with  expectations  of  local  accuracy  low. 

Horizontal  resolution  and  local  geometric  detail 
are  driven  by  two  concerns:  minimizing  deviations 


from  grid  orthogonality,  a  constraint  imposed  by 
the  ELCIRC  algorithm  (Zhang  et  al.,  2004);  and 
keeping  numerical  diffusion  under  control.  As 
discussed  in  detail  elsewhere  (Baptista  et  al.,  2004), 
we  have  already  achieved  orthogonality  for  most 
elements  in  the  grid  (Fig.  7b, d),  but  further  grid 
refinement  is  useful  to  reduce  numerical  diffusion 
(Fig.  11),  and  thus  to  enhance  the  ability  of 
ELCIRC  to  represent  coastal  eddies  and  other 
sharp  spatial  features. 

The  vertical  grid  consists  of  62  z-levels,  with 
finer  resolution  concentrated  on  the  top  30  meters 
of  the  water  column  (Fig.  8).  As  typical  in 
z-coordinate  models,  near-bottom  representation 
is  challenging  for  ELCIRC.  With  our  choice  of 
vertical  and  horizontal  grids,  difficulties  in  the 
estuary  and  near  plume  are  minimized,  at  the 
expense  of  the  continental  shelf,  continental  slope, 
and  deep  ocean  (Fig.  9). 

Capitalizing  on  the  ability  of  ELCIRC  to  handle 
Courant  numbers  well  above  unity,  a  time  step  of 
1.5  min  is  used.  With  this  time  step,  Courant 
numbers  as  large  as  4  (in  the  estuary;  Fig.  10)  and 
even  10  (upriver;  figure  not  shown)  are  common.  A 
time  step  of  as  large  as  15  min  would  still  have 
been  appropriate  for  purely  barotropic  simula¬ 
tions.  However,  time  steps  larger  than  1.5  min  lead 
to  parasitic  oscillations  in  the  vicinity  of  strong 
baroclinic  forcings. 

To  understand  these  oscillations,  we  note  that  a 
Courant-Friedrich-Levy  condition  associated 
with  the  baroclinic  term  can  be  estimated  via  the 
maximum  internal  wave  speed  as 


A  t^/cfh 
u  Ax 


<1, 


(2) 


where  g'  =  gAp/p0  is  the  reduced  gravity.  There¬ 
fore,  a  theoretical  maximum  time  step  for  stability 
can  be  estimated  as 


At  TV 
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S  79  s, 


(3) 


assuming  a  typical  channel  depth  of  h  =  20  m, 
a  horizontal  resolution  of  175  m  (e.g.,  see 
Fig.  6),  and  an  extreme  case  of  freshwater 
(p  =  1000  kg/m3)  at  the  surface  and  saltwater 
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Fig.  6.  Equivalent  diameters  near  the  Columbia  River  entrance  for  a  typical  CORIE  grid,  (a)  Isolines;  both  color  and  “relief’  indicate 
the  magnitude  of  the  equivalent  diameter,  (b)  Histogram,  (c)  A  zoom  into  the  histogram  for  equivalent  diameters  up  to  and  including 
2000  m. 


(p  =  1025  kg/m3)  at  the  bottom.  This  order-of- 
magnitude  estimate  lends  credence  to  the  trial  and 
error  analysis  that  led  to  the  choice  of  an 
operational  time  step  of  1.5  min. 

Based  on  the  formal  analysis  of  Casulli  and 
Cattani  (1994),  we  set  the  implicitness  factor  9 
(see  Eq.  (37)  in  Zhang  and  Baptista,  2004)  to  0.6 
for  best  accuracy,  thus  weighing  the  present 
time  step  slightly  more  than  the  previous  time 
step  in  treating  terms  of  the  continuity 
and  momentum  equations  that  are  handled 


implicitly.  Empirical  trial-and-error  supported  this 
choice. 

By  default,  in  ELCIRC,  the  tracking  of  char¬ 
acteristic  lines  is  performed  with  a  simple  Euler 
integration.  N  integration  steps  (referred  to  as  sub¬ 
time  steps)  are  allowed  for  tracking  between  times 
n  +  1  and  n.  In  DB 1 1 ,  we  let  the  sub-time  step  be 
chosen  automatically,  to  account  for  local  gradi¬ 
ents  of  velocity.  Specifically: 

N  =  max{Wm,n,  min  [Nmax,  ma\(Nx,  Ny,  N:)])  (4) 
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(c)  (d) 


Fig.  7.  Orthogonality  indices  0.  (a)  and  (c):  Definition  sketches  for  a  triangle  and  quadrangle,  (b)  and  (d):  Indices  for  triangles  and 
quadrangles.  For  orthogonal  triangles,  8>  0.  For  orthogonal  quadrangles,  9=1. 


where  Nmi„  and  Nmax  are  user-specified  limits  (2 
and  9,  for  DB11)  and 
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where  u,  v,  and  w  are  velocities  in  the  x,  y,  and  z 
directions;  horizontal  gradients  |6u/6/|  and  |0a/9/| 


are  computed  along  sides  of  elements;  and  A z 
is  the  local  vertical  grid  size.  Experiments  where 
we  enforced  smaller  sub-time  steps  or  tracked 
the  characteristic  lines  using  higher-order 
methods  (e.g.,  5th  order  Runge-Kutta)  revealed 
no  substantial  gains  in  accuracy  (results  not 
shown). 

Three  types  of  physical  parameterization  play 
potentially  important  roles  in  the  model  output: 
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Fig.  8.  Vertical  discreatization,  as  represented  by  Az,  is  shown  as  a  function  of  depth  (a),  including  a  detail  inset  for  the  top  30  m  (b). 


bottom  friction,  surface  stress,  and  vertical  mixing. 
Choices  for  DB11  were  as  follows: 

Bottom  friction'.  The  existence  of  different  bed 
forms  in  the  Columbia  River,  upstream  and 
downstream  of  the  Astoria-Tongue  Point  region, 
has  long  been  recognized  (Hamilton,  1990).  We 
coarsely  represent  this  difference  by  imposing  a 
spatially  varying  bottom  drag  coefficient  (Cob\  see 
Eq.  (14)  in  Zhang  et  al.,  2004).  Specifically,  we 
allow  for  a  frictionless  bottom  (Cob  =  0)  in  the 
continental  shelf  and  in  the  Columbia  River  up  to 
20  km  upstream  of  the  estuary  entrance,  we  impose 
substantial  friction  (Cm  =  0.0045)  above  30  km 
upstream  of  the  estuary  entrance,  and  we  let  Cm 
transition  linearly  between  these  two  regions. 
While  characterization  of  bottom  friction  is  not  a 
closed  issue  for  the  Columbia  River  (e.g.,  Frolov 
et  al.,  2004),  improvements  of  model  results  based 
on  optimizing  values  of  Cob  have  been  modest. 
Internal  calculations  of  Cm  based  on  matching 
model  velocities  with  bottom  boundary  layer 
profiles  (Eq.  (15)  in  Zhang  et  al.,  2004)  have  not 
proved  clearly  superior,  either,  possibly  because  of 
the  difficulty  of  z-coordinate  models  such  as 


ELCIRC  in  representing  the  bottom  boundary 
layer  (Fig.  11). 

Surface  stress :  Through  sensitivity  analysis, 
partially  reported  in  Section  3.2,  we  found  that 
the  bulk  aerodynamic  algorithm  of  Zeng  et  al. 
(1998)  to  be  superior  to  more  traditional,  simpler, 
but  less  process-driven  surface  stress  formulations 
(e.g.,  see  review  in  Pond  and  Pickard,  1998).  As 
described  in  Zhang  et  al.  (2004),  the  algorithm  of 
Zeng  et  al.  (1998)  accounts  for  surface  layer 
stability,  free  convection,  and  variable  roughness 
length  at  the  ocean-atmosphere  interface.  This 
algorithm  is  now  standard  in  the  CORIE  modeling 
system  whenever  simulations  (such  as  in  DB11) 
use  the  output  of  atmospheric  models  to  drive 
wind  fields  and  heat  balance  budgets. 

Vertical  mixing :  ELCIRC  offers  various  alter¬ 
natives  to  characterize  vertical  mixing  (see  Section 
2.4  of  Zhang  et  al.,  2004).  DB11  uses  a  k-kl 
closure,  with  mixing  limits  imposed  as  shown  in 
Table  5.  The  choice  of  different  mixing  limits  was 
meant  to  reflect  different  mixing  regimes  inside 
and  outside  the  estuary  and,  more  specifically,  to 
prevent  over-mixing  inside  the  estuary.  The  choice 
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Fig.  9.  Bottom  representation  in  the  numerical  model,  (a)  The  actual  bathymetry  in  ‘gray’  and  its  representation  in  the  model  in  ‘black’ 
along  a  transect  extending  from  the  estuary  to  offshore.  The  orientation  of  the  transect  is  shown  in  the  inset  map.  Distance  (km)  is  from 
the  origin  of  the  transect  inside  the  estuary,  (b)  The  difference  between  the  model  representation  and  actual  bathymetry.  Positive  values 
indicate  regions  where  the  model  bathymetry  is  shallower  than  actual  bathymetry.  (c)-(d)  The  same  transect  but  zoomed  inside  the 
estuary. 
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Fig.  10.  Estimated  Courant  number  near  the  Columbia  River  entrance  for  a  high-discharge  condition.  Both  color  and  “relief’  indicate 
the  magnitude  of  the  Courant  number,  which  approaches  4  in  the  region  shown. 
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Fig.  11.  Maximum  numerical  diffusivity  (m2/s)  near  the  Columbia  River  entrance.  Both  color  and  “relief’  indicate  the  magnitude  of 
the  numerical  diffusivity.  The  2D  inset  shows  a  larger  area,  generally  with  larger  numerical  diffusivity  because  of  lower  resolution. 
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has  significant  implications  for  the  results 
(see  Section  3.2). 

Arguably,  external  forcings,  in  particular  the 
ocean  conditions  and  atmospheric  forcings, 
are  currently  the  single  most  significant  source 
of  uncertainty  for  CORIE  simulations  (see 
Section  2.1  for  the  details  for  initial  and  boundary 
conditions).  In  DB11,  S,  T  approximately  300  km 
or  more  from  the  mouth  of  the  estuary  are  nudged 
to  NCOM  results,  with  a  maximum  nudging  factor 
of  5%  (see  Section  2.3). 

3.1.2.  Representation  of  water  levels 

We  first  consider  water  levels  at  a  single  station, 
Tongue  Point  (or  tpoin,  using  the  CORIE  5-digit 
terminology  for  field  stations).  A  station  of  the 
NOAA  CO-OPS  tidal  network  (Fig.  5),  Tongue 
Point  is  located  inside  the  estuary,  about  30  km 


upstream  from  the  entrance  (Fig.  4).  Long-term 
water-level  observations,  with  accurate  vertical 
data,  are  available;  the  record  for  2002,  in 
particular,  is  uninterrupted  (Fig.  12e). 

DB11  provides  a  robust  description  of  Tongue 
Point  water  levels.  The  average  error1  for  a  full- 
year  simulation  is  —0.01m,  with  a  standard 
deviation  of  0.188  m  and  a  root  mean  square  error 
of  0.188  m.  Water  levels  are  over-estimated  at 
most  by  0.658  m  and  are  underestimated  at 
most  by  0.848  m.  Histograms  of  errors  are  shown 
in  Fig.  13,  both  for  the  full  signal  and  for  specific 
frequency  bands,  defined  as  low  pass  (r>30h), 
band  pass  (9.6h<T ^30h;  i.e.,  in  the  astronomic 


'Throughout  this  paper,  model  error  is  defined  as  “simula¬ 
tions  minus  observations”.  This  definition  ignores  observation 
errors,  but  is  sufficient  for  the  purposes  of  our  discussion. 
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0.00 

0.00 

-0.01 
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0.03 

0.12 
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0.28 

0.05 

0.20 

0.22 
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-0.34 

-0.05 

-0.22 
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Maximum  over  prediction 

0.66 

0.13 

0.38 

0.30 

Maximum  under  prediction 

-0.85 

-0.14 

-0.52 

-0.75 

Fig.  13.  DB11  error  statistics  for  water  levels  at  Tongue  Point  (in  meters)  for  baseline  simulation,  (a)  Full  signal,  (b)  High-pass,  (c) 
Band-pass,  (d)  Low-pass. 


tidalrange),  and  high  pass  ( T  <9.6h;  i.e.,  inclusive 
of  shallow  water  tides). 

Time  series  of  errors  at  Tongue  Point 
(Fig.  12  b,  d)  suggest  that  band-pass  and  high- 
pass  errors  respond  directly  to  tidal  forcing  and 
tend  to  be  largest  during  spring  tides.  Band-pass 
errors  are  substantially  larger  than  high-pass 
errors,  reflecting  the  relative  difference  of  the 
signals  represented  (astronomic  tides  being  larger 
than  shallow  water  tides).  While  tides  are  non¬ 
stationary  in  the  Columbia  River  due  to  interac¬ 
tions  with  river  discharge,  harmonic  analysis  for 
the  full  year  of  2002  (Tables  1  and  2)  provides  a 
useful,  if  simplified,  context.  We  note,  in  particular 
that  in  most  constituents,  the  tidal  amplitudes  are 
over-estimated  and  the  model  simulations  lead  the 
data  (smaller  phase  lag).  The  exception  is  the  M6 
component,  where  the  amplitudes  are  under¬ 
predicted  by  the  model,  thus  suggesting  insufficient 
bottom  friction. 

Low-pass  errors  (Fig.  12c)  show  a  seasonal 
trend,  with  the  model  tending  to  overestimate 


observations  in  summer  and  to  underestimate 
them  in  winter.  Strong  winter  storms  in  January 
and  December  introduce  the  largest  errors  (see 
also  Section  3.2).  We  note  that  the  model  is  able  to 
generate  internally  a  significant  part  of  the  low- 
pass  signal,  even  though  that  signal  is  not  forced  at 
the  domain  boundaries.  Coastal  winds  and  atmo¬ 
spheric  pressure  gradients  are  responsible  for  this 
internal  generation. 

A  comparison  of  error  patterns  across  selected 
stations  of  the  NOAA  CO-OPS  tidal  network  is 
shown  in  Fig.  14  (see  Fig.  5  for  station  locations). 
There  is  remarkable  spatial  coherence  at  all 
frequencies,  with  some  informative  exceptions.  In 
particular: 

•  The  average  error  at  cnbw\,  a  station  at  the 
entrance  of  the  Strait  of  Juan  de  Fuca,  is 
substantially  larger  than  that  of  any  other 
station  (cf.  Table  3).  This  may  be  due  to  the 
overly  simplistic  representation  of  that  part  of 
the  domain  in  the  CORIE  modeling  system.  The 
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Table  1 

Illustrative  tidal  constituents  at  the  ocean  boundary  and  in  the  Columbia  River  estuary 


Frequency  (rad  s  ') 

Amplitude  imposed  at  ocean  boundary  (m)“ 

Amplitude  at  Tongue  Pointb  (m) 

Maximum 

Minimum 

46  °  1 5'N 

Data 

DB06 

m2 

1.405189e  —  04 

1.0064 

0.3395 

0.8518 

0.9255 

0.9987 

S2 

1.454441e  —  04 

0.2923 

0.0626 

0.2379 

0.2400 

0.2844 

n2 

1.378797e  —  04 

0.2063 

0.0815 

0.1733 

0.1798 

0.1963 

k2 

1.458423e  —  04 

0.0145 

0.0645 

0.1030 

0.1129 

K, 

7.2921 17e- 05 

0.4666 

0.3289 

0.4247 

0.4651 

0.5148 

0, 

6.759775e  -  05 

0.191 

0.2587 

0.3063 

0.3173 

Pi 

7.25 1 056e  —  05 

0.1461 

0.1030 

0.1317 

0.1158 

0.1326 

Q, 

6.495457e  -  05 

0.0353 

0.0455 

0.0504 

0.0556 

m4 

2.810378e  -  04 

— 

— 

— 

m6 

4.21 5567e  —  04 

— 

— 

— 

“The  ocean  boundary  of  the  CORIE  modeling  domain  extends  from  35  °N  to  50  °N.  46°15'N  is  the  latitude  of  the  Columbia  River 
entrance.  Values  are  extracted  from  Myers  and  Baptista  (2001). 

bTides  are  non-stationary  in  the  Columbia  River  because  of  the  strong  nonlinear  influences  of  river  discharge.  Amplitudes  shown  for 
Tongue  Point,  a  NOAA  tidal  station  approximately  30  km  upstream  of  the  mouth  of  the  Columbia  River,  are  from  harmonic  analysis 
for  year  2002.  DB06  refers  to  a  specific  CORIE  simulation  (see  text,  Section  3.1). 


Table  2 

Comparison  of  observed  and  simulated  tidal  constituents  at 
Tongue  Point  for  2002 


Amplitude  at  Tongue 
Point  (m) 

Phase  at  Tongue  Point 
(deg) 

Data 

DB1 1 

Data 

DB11 

m2 

0.9255 

0.9930 

159.60 

155.34 

s2 

0.2400 

0.2832 

293.85 

290.27 

n2 

0.1798 

0.1952 

336.53 

323.38 

k2 

0.1030 

0.1141 

310.11 

307.12 

K, 

0.4651 

0.5076 

267.87 

263.34 

0, 

0.3063 

0.3183 

117.13 

110.77 

Pi 

0.1158 

0.1316 

256.58 

249.91 

Q. 

0.0504 

0.0595 

316.51 

308.16 

m4 

0.0057 

0.0163 

249.72 

279.20 

m6 

0.0109 

0.0028 

173.10 

174.03 

effect  is  seen  clearly  in  the  low-pass  and  band¬ 
pass  errors  (Fig.  14a,  c)  and  in  the  amplitudes  of 
specific  tidal  constituents  (in  particular  for  M2 
and  M4;  Table  4). 

•  Low-pass  errors  (Fig.  14a)  tend  to  show  near- 
synchronous  spikes  (underestimations)  during 
strong  winter  storms,  consistent  with  the  regio¬ 
nal  development  of  frontal  systems.  The  south¬ 


ernmost  station  (cmoc  1)  is  the  least  correlated  of 
all  stations.  A  detailed  analysis  of  the  correlation 
between  frontal  systems  and  water-level  re¬ 
sponses  will  be  presented  in  a  separate  article. 

•  High-pass  errors  (Fig.  14b)  are  typically  small 
and  smaller  than  band-pass  errors.  The  two 
cases  {cwbwl  and  skaw\)  where  substantial  high- 
pass  errors  occur  are  in  stations  located  in 
shallow  areas  with  very  poor  grid  and  bathy¬ 
metric  resolution  (in  Willapa  Bay  and  in  a 
freshwater  confluence  of  the  Columbia  River, 
respectively).  High-pass  errors  at  these  stations 
are  dominated  by  shallow  water  tidal  frequen¬ 
cies,  an  indication  that  local  resolution  is 
insufficient  to  account  for  strong  tidal  nonlinea¬ 
rities  occurring  in  these  areas. 

•  Errors  at  Tongue  Point  are  generally  coherent 
with  dominant  errors  at  coastal  stations.  How¬ 
ever,  the  spring-neap  asymmetry  of  band-pass 
errors  at  Tongue  Point  is  larger  than  at  coastal 
stations  (except  for  (cnbw  1  and  cwbwl),  which 
are  affected  by  special  circumstances). 

3.1.3.  Representation  of  wetting  and  drying 
The  Columbia  River  is  subject  to  extensive 

wetting  and  drying.  Representation  of  this  process 
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is  challenging,  but  necessary  and  within  the  type  of 
capabilities  expected  of  ELCIRC.  In  Fig.  15,  we 
provide  a  snapshot,  during  low-water  conditions, 
of  a  SAR  image  that  identifies  (with  some 
subjectivity)  wet  and  dry  areas,  and  a  DB11 
simulation  displaying  equivalent  information. 
For  reference,  we  also  include  a  representation 
of  the  DB11  numerical  grid,  which  permits 
identification  of  the  areas  that  are  kept  perma¬ 
nently  dry  in  the  modeling  domain.  The  three 
images  (Fig.  15a-c)  are  consistently  geo-refer¬ 
enced,  and  the  SAR  and  DB1 1  images  are  synoptic 
within  5  min. 

The  qualitative  agreement  between  the  simula¬ 
tions  and  the  SAR  image  is  remarkably  good  both 


in  the  mainstem  of  the  estuary  and  in  the 
Cathlamet  Bay.  We  found  the  results  very 
encouraging  relative  to  the  ability  of  ELCIRC  to 
represent  wetting  and  drying  processes. 

3.1.4.  Representation  of  salinity  fields 
Columbia  River  salinity  fields  vary  dramatically 
in  space  and  time  (Jay  and  Smith,  1990),  with 
major  oceanographic  and  ecological  implications. 
Multi-scale  representation  of  this  variability  is 
extremely  challenging  and  is  a  major  goal  of  the 
CORIE  modeling  system.  This  paper  provides 
only  an  introduction  to  challenges  and  successes 
towards  this  goal.  We  begin  by  examining  year¬ 
long  time  series  of  salinities  at  three  stations, 


Fig.  14.  DB1 1  errors  in  water  levels  at  multiple  tide  gauges  along  the  coast  and  in  the  Columbia  River  system  (see  Figs.  4  and  5  for 
locations):  (a)  Low-pass,  (b)  High-pass,  (c)  Band-pass. 
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Fig.  14.  ( Continued ) 


where  CORIE  sensors  are  available:  ogiOl,  red26, 
and  am\(>9. 

Station  ogiOl  is  located  in  the  continental  shelf, 
at  a  depth  of  about  100  m,  and  is  approximately 
25  km  southwest  of  the  mouth  of  the  Columbia 
River  estuary  (Fig.  4).  The  station  is  in  the  path  of 
the  near  plume  during  periods  of  northerly  wind, 
which  are  dominant  during  the  spring  and 
summer.  During  the  winter,  southerly  winds 
dominate,  and  ogiOl  is  off  the  plume  path  for 
extended  periods  (e.g.,  see  pattern  of  observed 
daily  maximum  and  minimum  salinities  in  January 
and  February,  Fig.  16a).  DB11  simulations 


(Fig.  16b)  show  a  plume  that  is  responsive  to 
changes  in  wind  direction. 

Station  red26  (Fig.  4)  is  approximately  13  km 
upstream  of  the  mouth  of  the  estuary,  near  the 
navigation  channel.  Bottom  daily  maximum  sali¬ 
nities  from  DB11  track  well  those  salinities 
observed  at  the  bottom  CORIE  sensor  (Fig.  17). 
In  particular,  we  note  the  good  correlation 
between  observed  and  simulated  sudden  dips  in 
salinity  (e.g.,  at  the  end  of  January  and  beginning 
of  July).  The  conditions  that  generate  these  dips 
are  under  investigation,  one  hypothesis  being  that 
the  dips  coincide  with  the  simultaneous  occurrence 
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Fig.  14.  ( Continued ) 


Table  3 

Comparison  of  error  statistics  for  water  levels  at  multiple 
NOAA  CO-OPS  for  DB11 


Stations 

Error  (m) 

Average 

95%  percentiles 

Std.  Dev. 

RMS 

cnbwl 

0.36 

0.03 

0.67 

0.20 

0.41 

cwbw\ 

0.04 

-0.32 

0.43 

0.23 

0.24 

tpoin 

-0.01 

-0.34 

0.28 

0.19 

0.19 

jteu’l 

0.07 

-0.27 

0.34 

0.19 

0.20 

csboi 

-0.05 

-0.35 

0.23 

0.18 

0.18 

cccc  1 

0.10 

-0.11 

0.30 

0.13 

0.16 

cnsc  1 

-0.03 

-0.27 

0.21 

0.15 

0.16 

caccl 

-0.03 

-0.24 

0.17 

0.14 

0.14 

cmoel 

-0.05 

-0.22 

0.11 

0.10 

0.11 

of  downwelling-favorable  winds  and  neap  tides. 
Like  those  observed,  daily  minimum  salinities 
from  DB11  decrease  substantially  during  spring 
freshets.  In  general,  however,  a  comparison  of 
observed  and  simulated  salinities,  and  daily  mini¬ 
mum  salinities  suggests  that  the  simulated  estuary 
is  fresher  during  ebb  than  it  should  be,  particularly 
during  the  latter  half  of  the  year. 

Station  am  169  (Figs.  4  and  18)  is  located  further 
upstream  in  the  estuary,  approximately  20  km 
from  the  mouth  and  also  near  the  navigation 
channel.  Often  fresh  during  ebb,  and  occasionally 
fresh  throughout  the  entire  tidal  cycle  during  spans 
of  the  spring  freshet,  amb\69  is  a  challenging 
station  to  represent  numerically.  In  particular, 
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Table  4 

Tidal  amplitudes  (m)  at  selected  NOAA  CO-OPS  stations 


caccl 

ccc\ 

cc/io3 

cmocl 

Data 

DB1 1 

Data 

DB1 1 

Data 

DB1 1 

Data 

DB1 1 

m2 

0.5686 

0.5186 

0.7103 

0.7003 

0.8060 

■SB 

0.4427 

S2 

0.1382 

0.1259 

0.1825 

0.1798 

0.2151 

0.1222 

n2 

0.1240 

0.1082 

0.1522 

0.1430 

0.1104 

0.1028 

k2 

0.0419 

0.0337 

0.0560 

0.0508 

0.0623 

0.0610 

0.0409 

0.0322 

K, 

0.3859 

0.4039 

0.4068 

0.4356 

0.4145 

0.4526 

0.3803 

0.3921 

0, 

0.2467 

0.2338 

0.2560 

0.2548 

0.2591 

0.2681 

0.2430 

0.2234 

P. 

0.1167 

0.1183 

0.1238 

0.1275 

0.1257 

0.1331 

0.1157 

0.1135 

Qi 

0.0438 

0.0446 

0.0457 

0.0489 

0.0468 

0.0496 

0.0432 

0.0419 

m4 

0.0006 

0.0002 

0.0007 

0.0007 

0.0097 

0.0003 

0.0016 

0.0001 

m6 

0.0005 

0.0001 

0.0004 

0.0003 

0.0040 

0.0004 

0.0002 

0.0001 

cnbw\ 

cnsc\ 

cpoo7> 

csboh 

Data 

DB11 

Data 

DB1 1 

Data 

DB1 1 

Data 

DB11 

m2 

0.7807 

1.0755 

0.6815 

0.6488 

0.7452 

0.7492 

0.8837 

S2 

0.2309 

0.3125 

0.1722 

0.1624 

0.1963 

0.1978 

0.2422 

n2 

0.1672 

0.2196 

0.1455 

0.1338 

0.1577 

0.1525 

0.1870 

k2 

0.0631 

0.0907 

0.0497 

0.0453 

0.0594 

0.0578 

0.0715 

K, 

0.5126 

0.5028 

0.4088 

0.4289 

0.4432 

0.4495 

0.4521 

0.4662 

O, 

0.3210 

0.3056 

0.2579 

0.2504 

0.2824 

0.2633 

0.2798 

0.2778 

p. 

0.1564 

0.1461 

0.1247 

0.1270 

0.1352 

0.1270 

0.1367 

0.1375 

Q. 

0.0574 

0.0532 

0.0456 

0.0475 

0.0498 

0.0504 

0.0503 

0.0502 

m4 

0.0104 

0.0013 

0.0117 

0.0011 

0.0008 

0.0002 

0.0136 

0.0010 

m6 

0.0073 

0.0002 

0.0062 

0.0007 

0.0013 

0.0001 

0.0072 

0.0002 

cwbw  1 

Data 

DB11 

m2 

0.9612 

0.8174 

S2 

0.2583 

0.2249 

n2 

0.1970 

0.1636 

k2 

0.0780 

0.0761 

K, 

0.4449 

0.4481 

0, 

0.2749 

0.2602 

Pi 

0.1373 

0.1288 

Qi 

0.0475 

0.0456 

m4 

0.0139 

0.0465 

m6 

0.0083 

0.0046 

DB11  shows  consistently  lower  levels  of  salt  than 
those  observed,  including  a  tendency,  for  most  of 
the  year,  to  predict  freshwater  conditions  during 
ebb,  even  at  the  bottom  of  the  estuary. 

The  results  presented  above  (and  additional 
results  available  at  CCALMR,  1996-2004)  suggest 
that,  on  the  whole,  DB11  salinities  are  responsive 


to  river  forcings,  coastal  winds,  and  spring-neap 
cycles,  but  that  salt  penetration  in  the  estuary 
tends  to  be  under  predicted.  As  a  baseline  for 
further  analysis  in  Section  3.2,  year-long  salt 
volumes  in  the  estuary  (in  the  form  of  15-day 
running  averages)  and  plume  volumes  (defined  by 
the  30psu  contour)  are  plotted  in  Fig.  19  against 
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Fig.  16.  (a)  Observed  daily  maximum  (red)  and  minimum  (blue)  salinity  at  the  CORIE  station  ogiOl,  from  a  sensor  located  1  m  below 
the  water  surface,  (b)  Daily  maximum  (red)  and  daily  minimum  (blue)  salinity  at  the  surface  from  DB11.  (c)  Vertical  structure  of 
salinity  at  ogiOl,  top  30m  from  baseline  simulations,  (d)  Daily  tidal  range  from  observations  at  Tongue  Point,  (e)  N-S  component  of 
daily-averaged  wind  speed  from  observations  at  the  NOAA  buoy  46029.  (f)  Daily-averaged  river  discharge  measured  at  Bonneville 
Dam.  Time  is  denoted  by  MM/YY. 


the  river  discharge  at  Bonneville  Dam.  The  salt 
volume  in  the  estuary  behaves  qualitatively  as 
expected,  responding  to  both  the  seasonal  varia¬ 
tion  of  river  discharges  and  to  spring-neap  cycles. 
Unlike  earlier  CORIE  databases,  there  are  no 
discontinuities  among  simulation  ensembles  in 
DB11,  suggesting  that  dynamic  equilibrium  has 
been  reached  within  each  ensemble.  The  lack  of 


discontinuities  validates  our  strategy  of  temporal 
parallelization  and  suggests  that  the  “response 
time”  of  the  system  is  on  the  order  of  one  to  two 
months.  Other  sensitivity  runs  revealed  that  the 
more  realistic  turbulence  closure  scheme  of  k-kl 
and  the  nudging  to  the  NCOM  results  in  the  ocean 
are  responsible  for  the  quick  convergence  of  flow 
fields  (not  shown  here). 
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Fig.  17.  (a-c)  Observed  daily  maximum  (red)  and  daily  minimum  (blue)  salinity  at  the  CORIE  station  red26  from  sensors  located  at 
3.3,  7.5  and  9.0m  below  MSL,  respectively,  (d)  Daily  maximum  (red)  and  daily  minimum  (blue)  salinity  at  red26,  bottom  layer,  from 
DB1 1 .  (e)  Vertical  structure  of  salinity  at  red 26  from  baseline  simulations,  (f)  Daily  tidal  range  from  observations  at  Tongue  Point,  (g) 
N-S  component  of  daily-averaged  wind  speed  from  observations  at  the  NOAA  46029  buoy,  (h)  Daily-averaged  river  discharge 
measured  at  Bonneville  Dam.  Time  is  denoted  by  MM/YY. 


Besides  fixed  stations,  scientific  cruises  are  also 
very  useful  in  gauging  the  model  performance  by 
covering  a  large  area.  A  comprehensive  survey  was 
conducted  in  May-July  2004,  which  generated  a 
wealth  of  observational  data.  Shown  in  Figs.  20 
and  21  are  some  preliminary  data-model  (from 
CORIE  daily  forecasts)  comparisons:  the  good 


qualitative  agreement  between  observed  and  simu¬ 
lated  plume  in  Fig.  21  shows  how  dynamic  the 
plume  is  to  the  forcing  within  several  days,  while 
Fig.  20  suggests  that  the  model  qualitatively 
captures  some  key  characteristics  of  the  plume, 
although  the  modeled  plume  tends  to  be  fresher 
than  the  observed  one,  especially  near  the  surface. 
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Fig.  18.  (a-c)  Observed  daily  maximum  (red)  and  daily  minimum  (blue)  salinity  at  the  CORIE  station  am  1 69  at  sensors  located  at 
2.6m,  11.3  m,  and  14.3  m  below  MSL,  respectively,  (d)  Daily  maximum  (red)  and  daily  minimum  (blue)  salinity  at  am\69,  bottom 
layer,  from  DB11.  (e)  Vertical  structure  of  salinity  at  am\69  from  baseline  simulations.  (0  Daily  tidal  range  from  observations  at 
Tongue  Point,  (g)  N-S  component  of  daily-averaged  wind  speed  from  observations  at  the  NOAA  buoy  46029.  (h)  Daily-averaged  river 
discharge  measured  at  Bonneville  Dam.  Time  is  denoted  by  MM/YY. 


A  detailed  report  on  the  May-July  2004  surveys 
will  be  published  elsewhere. 

3.2.  Sensitivity  to  modeling  choices 

To  explore  the  sensitivity  of  the  results  to 
modeling  choices,  we  present  in  this  section 


selected  results  from  various  databases  and  cali¬ 
bration  runs  as  defined  in  Table  5. 

3.2.1.  Parameterization  of  surface  stresses 

Figure  22  motivated  our  preference  for  a 
particular  parameterization  of  surface  stresses 
(Zeng  et  al.,  1998).  This  figure  shows  the 
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Fig.  19.  Integral  salt  balance  metrics  computed  from  DB1 1.  (a)  Volume  of  salt  inside  the  estuary,  (b)  Plume  volume  based  on  the  30psu 
contour  as  defined  for  plume  extent.  Note  that  there  is  no  visible  sharp  discontinuities  across  ensembles,  (c)  Daily-averaged  river 
discharge  measured  at  Bonneville  dam. 


week-averaged,  root-mean-square  errors  of  low- 
pass  water  levels  at  Tongue  Point,  contrasted 
against  the  x-component  (i.e.,  approximately,  the 
E-W  component)  of  the  Ekman  transport  in  the 
vicinity  of  the  Columbia  River.  The  x-component 
of  the  Ekman  transport  is  computed  from  external 
forcings  and  ELCIRC  results  as  (following  Cush- 
man-Roisin,  1994): 


E%  =  J  dt  J  U(x,y)ds  (in  m3), 


U{x,y) 


i: 


u(x,y,z)dz 


*y(x,y) 

Pof 


(in  m2/s), 


(6) 


where  L  bounds  a  rectangle  centered  at  the  mouth 
of  the  estuary,  T  =  7  days,  u  is  the  x-component  of 
the  water  velocity,  iy  is  the  ^-component  of  the 
wind  stress,  and  p0  and  /  are  respectively  a 
reference  density  and  the  Coriolis  factor. 

DB04  and  DB05  were  built  with  an  early  version 
of  ELCIRC  (version  4.01,  rather  than  5.01  used  in 
most  recent  runs).  Like  DB06,  they  were  initialized 
from  a  Levitus  condition  and  used  a  zero-equation 
parameterization  of  vertical  mixing  (Pacanowski 
and  Philander,  1981;  hereafter,  “P&P”),  although 
other  set-up  details  differ  from  DB06.  Most 
importantly,  DB05  is  the  first  hindcast  database 
to  use  the  approach  of  Zeng  et  al.  (1998)  to 
parameterize  surface  stresses,  while  DB04  uses  one 
of  the  empirical  relationships  of  Pond  and  Pickard 
(1998). 
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Fig.  20.  Comparison  of  model  results  to  data  for  an  offshore  cruise  conducted  on  July  15,  2004.  Thermosalinigraph  (TSG)  data  were 
collected  along  the  entire  path  (a),  and  conductivity-salinity-depth  (CTD)  casts  were  performed  at  multiple  locations  (marked  by  +). 
Model  and  observed  data  were  referenced  relative  to  the  free  surface,  (b)  Model  data  (DB11)  from  1  m  below  the  free  surface  (blue) 
compared  to  TSG  data  (red),  (c)  Salinity  isolines  constructed  from  CTD  cast  data  (location  varies  over  time),  (d)  Salinity  isolines  from 
DB11  along  cruise  path. 


DB04  responds  to  downwelling  events  early  and 
late  in  the  year  with  the  type  of  large  root  mean 
square  errors  that  were  characteristic  of  early 
databases  and  calibration  runs,  all  of  which  used 
Pond  and  Pickard  (1998).  In  DB05,  06,  and  11, 
this  type  of  errors  is  greatly  attenuated  relative  to 
DB04.  None  of  the  substantial  modifications 
introduced  after  DB05  (in  code  formulation, 
external  forcings,  or  simulation  set-up)  signifi¬ 
cantly  affects  the  response  of  errors  to  down- 


welling  conditions,  thus  confirming  surface  stress 
parameterization  as  the  transformative  element  in 
improving  in-estuary  responses  to  coastal  winds 
during  downwelling. 

The  remaining  errors  in  water  levels  during 
downwelling  regimes  are  currently  being  investi¬ 
gated.  The  prevailing  hypothesis  is  that  errors 
derive  predominantly  from  the  difficulty  of  ex¬ 
ternal  forecasts  to  represent  the  intensity  and  the 
phase  of  strong  winter  frontal  systems. 
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ctsy:  Dr.  Wesson,  N 


Fig.  21.  Comparison  of  plume  profiles  for  3  days  in  May  2004,  between  (a)  model  and  (b)  observation  from  an  airborne  survey 
(courtesy  of  Dr.  Wesson). 
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Table  5 

Summary  of  runs 


ELCIRC 

version 

Ocean  initial 
conditions 

Nudging  for 
ocean  S,T 

Wind  stress 
formulation 

Turbulence 
closure  scheme 

Diffusivity  limits  in 
the  ocean  (for  k-kl 
only)  (m2/s) 

Diffusivity  limits  in 
the  estuary  (for  k-kl 
only) 

DB04 

4.01 

Levitus 

No 

Pond  and  Pickard 

Pacanowski 

N/A 

N/A 

(1998) 

and  Philander 

(P&P)  (1981) 

DB05 

4.01 

Levitus 

No 

Zeng  et  al.  (1998) 

P&P 

N/A 

N/A 

DB06 

5.01 

Levitus 

No 

Zeng  et  al.  (1998) 

P&P 

N/A 

N/A 

DB10 

5.01 

NCOM 

Yes 

Zeng  et  al.  (1998) 

k-kl 

umax  =  10, 

umax  = 

10-\ 

umin  =  10-6 

umin  = 

10"6 

DBl  1“ 

5.01 

NCOM 

Yes 

Zeng  et  al.  (1998) 

k-kl 

umax  =  10, 

umax  = 

10-3, 

umin  =  10-6 

umin  = 

10~6 

C88a 

5.01 

NCOM 

Yes 

Zeng  et  al.  (1998) 

P&P 

umax  =  10, 

umax  = 

10-\ 

umin  =  10~6 

umin  = 

10~6 

ClOOa 

5.01 

S  =  32psu, 

Yes 

Zeng  et  al.  (1998) 

k-kl 

umax  =  10, 

umax  =s 

10-2, 

T=  15°C 

umin  =  10~4 

umin  = 

10-6 

ClOOb 

5.01 

S  =  35psu, 

Yes 

Zeng  et  al.  (1998) 

k-kl 

umax  =  10, 

umax  = 

10“2, 

T=  15°C 

umin  =  10~4 

umin  = 

10~6 

ClOOc 

5.01 

S  =  35psu, 

Yes 

Zeng  et  al.  (1998) 

k-kl 

umax  =  10, 

umax  = 

0.0005, 

T=  15°C 

umin  =  10~4 

umin  — 

10~6 

ClOOe 

5.01 

S  =  35psu, 

Yes 

Zeng  et  al.  (1998) 

k-kl 

umax  =  10, 

umax  = 

10-2, 

T=  15°C 

umin  =  10~3 

umin  = 

10~6 

“The  main  differences  between  DB\0  and DBl  1  are  that  the  latter  includes  the  Strait  of  Geogia  and  the  Fraser  River,  and  uses  more 
recent  bathymetric  information. 


3.2.2.  Salinity  in  the  estuary  and  the  plume 

The  transport  of  salt  and  heat  is  very  sensitive  to 
vertical  mixing.  While  early  databases  and  calibra¬ 
tion  runs  used  P&P,  mainly  due  to  its  efficiency, 
this  zeroth  order  closure  did  not  give  satisfactory 
results  in  either  the  plume  or  the  estuary.  To 
demonstrate  this,  results  from  DB10  (which  uses 
k-kl)  and  C88a  (which  uses  P&P)  (cf.  Table  5)  are 
compared;  the  only  difference  between  them  is  the 
turbulence  closure  scheme  used.  The  results  of 
C88a  are  very  different  from  DB10,  in  terms  of 
plume  shape  and  extent  (Fig.  23),  and  salinity-time 
series  inside  the  estuary  (Fig.  24).  The  overly  fresh, 
large,  and  non-responsive  plume  shown  in  Fig.  23b 
is  typical  of  all  runs  using  P&P  closure. 

While  DB11  uses  k-kl,  it  still  under-predicts  the 
extent  of  salt  intrusion  in  the  estuary.  Hence, 
further  test  runs  were  done  to  improve  the  results. 
Of  all  choices  of  model  parameters,  we  were 
particularly  interested  in  two:  the  “background” 
ocean  salinity,  and  the  mixing  limits  inside  and 


outside  the  estuary.  Studying  the  sensitivity  of  the 
simulations  to  the  first  parameter  was  aimed  at 
understanding  the  impact  of  underpredicted 
NCOM  salinities  in  the  continental  shelf  region. 
Tests  focused  on  the  second  parameter  were 
designed  to  understand  the  sensitivity  of  the 
simulations  to  detailed  choices  within  the  k-kl 
closure  scheme.  To  simplify  the  task,  we  opted  for 
a  non-stratified  ocean  of  constant  S,  T.  All 
sensitivity  runs  (ClOOa-e)  were  cold-started  from 
week  27  of  2004  and  covered  3  consecutive  weeks. 
While  details  of  the  runs  are  in  Table  5,  we 
summarize  here  their  major  differences: 

1.  Cl  00b:  Base  run,  with  k-kl  closure,  the  mini¬ 
mum  diffusivity  in  the  ocean  being  10_4m2/s, 
the  maximum  diffusivity  in  the  estuary  being 
10_2m2/s,  and  an  ocean  salinity  of  35psu; 

2.  ClOOa:  ClOOb  with  an  ocean  salinity  of  32psu; 

3.  ClOOc:  ClOOb  with  the  maximum  diffusivity  in 
the  estuary  being  0.0005  m2/s; 
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Fig.  22.  RMS  errors  (computed  over  1-week  intervals)  in  water  levels  at  Tongue  point,  (a)  Errors  in  DB04.  (b)  Errors  in  DB05.  (c) 
Errors  in  DB06.  (d)  Errors  in  DB11.  (e)  Ekman  transport.  Time  is  in  weeks. 


4.  ClOOe:  ClOOb  with  the  minimum  diffusivity  in 

the  ocean  being  10~3m2/s. 

Comparisons  were  made  in  week  29  for  the 
bottom  salinity  at  the  station  aw  169  (Fig.  25),  as 
well  as  the  instantaneous  intrusion  length  (from 
the  mouth  of  the  estuary)  in  the  south  channel, 
using  5psu  as  the  cut-off  value  (Fig.  26).  Cl 00c, 
with  a  larger  ocean  salinity  to  start  with  and 
smaller  maximum  diffusivity  in  the  estuary, 
provides  the  best  results  in  all  regards,  and 
Cl 00a,  with  a  smaller  ocean  salinity,  provides  the 
worst  results.  This  finding  confirms  that  the  ocean 
supply  of  salt  is  important  for  the  salinity  intrusion 
in  the  estuary.  ClOOe  predicts  a  deeper  intrusion 
but  a  lower  bottom  salinity  than  ClOOb,  indicating 
that  ocean  “ambient”  mixing  also  influences  the 
estuary.  Larger  ambient  mixing  in  the  ocean  would 
entrain  more  salt  water  but  mix  it  more  effectively, 
thus  being  consistent  with  both  reduced  bottom 
salinity  and  stronger  salt  intrusion.  Results  suggest 
an  under-prediction  of  ocean  mixing  in  typical 


CORIE  simulations,  which  might  relate  to  under¬ 
representation  of  wind  stress,  or  to  the  neglect  of 
global-scale  ocean  circulation. 

Note  that  even  the  best  results  (ClOOe)  in  the 
above  experiment  show  underpredicted  salinity  at 
upstream  stations  like  aw  169.  The  same  trend  is 
generally  true  for  experiments,  not  reported  here, 
involving  moderate  refinements  of  horizontal  and 
vertical  grids  inside  the  estuary,  sensitivity  to 
bottom-drag  coefficients,  and  updated  bathymetry 
according  to  more  recent  surveys.  This  persistent 
trend  suggests  that  ELCIRC  might  be  approach¬ 
ing  a  natural  limit  in  its  ability  to  represent  salt 
propagation  without  data  assimilation. 

4.  Final  considerations 

CORIE  offers  an  early  example  of  the  inclusion 
of  sustained  multi-scale  modeling  in  ocean  ob¬ 
servatories.  Through  systematic  experimentation, 
we  have  made  substantial  progress  towards  a 
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and  (b)  C88a. 


A.M.  Baptista  et  id.  /  Continental  Shelf  Research  25  (2005)  935-972 


967 


.0  L- SJg. _ l.  _ 1  _ l _ ft  .  «  w  .  i  v  w  ^  h  l 

5-29  5-30  5-31  6-1  6-2  6-3  6-4  6-5 


Date 

Fig.  24.  Comparison  of  bottom  salinity  at  tansy  in  week  22  of  2002.  Note  that  there  is  a  gap  in  data. 


physically  based  description  of  the  baroclinic 
circulation  of  the  Columbia  River  estuary  and 
plume. 

ELCIRC  has  been  an  integral  part  of  that 
progress.  Indeed,  the  introduction  of  ELCIRC 
into  the  CORIE  modeling  system  in  2001  has  been 
transformative  of  our  ability  to  conduct  a  very 
large  number  of  meaningful,  long-term  3D  bar¬ 
oclinic  simulations.  Together  with  extensive  long¬ 
term  observations,  these  simulations  are  providing 
new  insights  into  the  Columbia  River  circulation. 

While  we  will  further  develop  the  theme  else¬ 
where,  the  modeling  work  reported  here  already 
shows  how  responsive  the  Columbia  River  is  to 
continental  shelf  processes — thus  opening  the 
doors,  for  instance,  to  understanding  the  impact 
of  climatic  cycles  on  the  estuary  and  plume,  such 
as  El  Nino-Southern  Oscillation  and  Pacific 
Decadal  Oscillation. 

From  our  perspective,  two  challenges  stand  in 
the  critical  path  of  CORIE  as  an  effective 
modeling  infrastructure  for  the  Columbia  River 
system.  One  challenge  relates  to  the  computational 
performance  of  ELCIRC.  Even  with  the  efficiency 


of  its  serial  version,  our  ultimate  goal  of  building 
multiple-decade  simulation  databases  will  require 
the  parallelization  of  ELCIRC,  an  on-going 
process.  The  availability  of  a  parallel  ELCIRC, 
coupled  with  the  40  CPUs  available  in  the  CORIE 
computer  cluster,  is  expected  to  significantly 
change  the  short-term  ability  of  the  CORIE 
modeling  system  to  generate  multi-year  hindcast 
databases. 

The  other  challenge — which  we  have  partially 
addressed  in  this  paper — relates  to  an  improved 
description  of  salt  dynamics  in  the  Columbia  River 
estuary  and  plume.  We  have  already  shown  that 
substantial  improvements  result  from  the  use  of  a 
2.5-equation  turbulence  closure  and  the  use  of 
ocean  conditions  derived  from  the  global  NCOM 
model.  However,  the  inability  of  ELCIRC  to 
resolve  the  bottom  boundary  layer  and  the  low- 
order  nature  of  the  ELCIRC  algorithm  are 
obstacles  that  need  to  be  addressed  moving 
forward. 

Two  options  are  being  considered  within  COR¬ 
IE:  one  involves  adding  data  assimilation  to 
ELCIRC,  while  the  other  involves  using  a  newly 
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(a)  Month-day  2004 


W  Month-day  2004 

Fig.  25.  Comparison  of  bottom  salinity  at  am\69  from  C100  a-c,  and  e.  (a)  ClOOb.  (b)  Differences  between  ClOOb  and  others.  Note 
that  the  true  bottom  salinity  should  be  much  higher  for  this  low  discharge  period. 


developed  higher-order,  terrain-following  coordi¬ 
nate  model.  While  the  application  of  the  new 
model  has  yielded  encouraging  preliminary  results 


for  the  salinity  intrusion  in  the  Columbia  River, 
the  eventual  path  forward  will  likely  be  determined 
by  computational  considerations.  In  any  case,  the 
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(a)  Month-day  2004 


(b)  Month-day  2004 


Fig.  26.  Comparison  of  salinity  intrusion  limits  based  on  5psu  threshold  for  (a)  ClOOb.  (b)  Differences  between  ClOOb  and  ClOOa,  c, 
and  e. 


robustness,  efficiency,  and  versatility  of  ELCIRC 
have  already  allowed  enormous  progress  in  under¬ 
standing  Columbia  River  circulation,  suggesting 
that  ELCIRC  is  a  valuable  tool  for  multi-scale 
circulation  modeling. 


5.  Appendix.  Definition  of  metrics  for  grid 
orthogonality 

Following  Casulli  and  Zanolli  (1998),  a  grid  is 
defined  as  orthogonal  if,  within  each  element,  a 
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point  (“center”,  although  not  necessarily  the 
geometric  center)  can  be  identified  such  that  the 
segment  joining  the  centers  of  two  adjacent 
elements,  and  the  side  shared  by  the  two  elements, 
have  a  non-empty  intersection  and  are  perpendi¬ 
cular  to  each  other. 

The  indices  defined  in  this  section  are  an  attempt 
to  provide  a  practical,  quantitative  metric  with 
which  to  evaluate  the  extent  to  which  hybrid  grids 
meet  the  orthogonality  requirement.  For  triangles, 
we  define  the  index  of  orthogonality  as 

(7) 

where  R  is  the  circum-radius  of  the  triangle,  and 
Lmjn  is  the  minimum  signed  distance  from  the 
circum-center  to  the  three  sides  (Fig.  7a).  The 
element  is  orthogonal  if  $3>0  (otherwise  non- 
orthogonal)  and  is  equilateral  if  S3  =  1 . 

For  quadrangles,  we  define  the  index  of 
orthogonality  as  (Fig.  7c): 

S4=^,  (0<94<oo) 

where  R  is  the  circum-radius  of  the  triangle  formed 
by  nodes  1-3, 2  and  J?4  is  the  distance  from  node  4 
to  the  circum-center  of  nodes  1-3.  The  element  is 
orthogonal  if  94  =  1.  Otherwise,  it  is  non-ortho- 
gonal.  Note  that  this  index  for  quadrangles 
assumes  that  the  circum-center  is  inside  the 
element,  a  requirement  that  first  needs  to  be 
checked. 
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